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Abstract 

Herein, we present analytical solutions for the electronic energy eigenvalues of the hydrogen molecular 
ion H 2, namely the one-electron two-fixed-center problem. These are given for the homonuclear case 
for the countable infinity of discrete states when the magnetic quantum number m is zero i.e. for 2 S + 
states. In this case, these solutions are the roots of a set of two coupled three-term recurrence relations. 
The eigensolutions are obtained from an application of experimental mathematics using Computer 
Algebra as its principal tool and are vindicated by numerical and algebraic demonstrations. Finally, 
the mathematical nature of the eigenenergies is identified. 



PACS: 31.15.-p, 31.15.Ar, 02.70. Wz, 31.50.Bc, 31-50.Df 



1 Introduction 



Although, there are well established software packages in the area of quantum chemistry such as GAUS- 
SIAN [1], MOLPRO [2] and GAMESS [3] which allow to obtain approximate numerical solutions to a 
number of fair sized molecules, the simplest molecule namely the hydrogen molecular ion, a quantum 
mechanical three-body problem, still remains mathematically intractable. 

In the fixed nuclei approximation, it is well known that the Schrodinger wave equation - a second order 
partial differential equation (PDE) - of the problem of one electron moving in the field of two fixed nuclei 
can be separated in prolate-spheroidal coordinates [4]. These coordinates allow a separation of variables 
that results in two non-trivial ordinary differential equations (ODE), and hence two eigenparameters: the 
energy parameter p 2 , and a separation constant A related to the total orbital angular momentum and the 
Runge-Lenz vector. 



We note that asymptotic expansions for small or large internuclear distances R have been obtained. A 
very comprehensive presentation of the energy eigenvalues for the ground state and a number of exited 
states is shown in the work of Cfzek et al. [9]. These could almost constitute analytical solutions but 
the resulting series are divergent though asymptotic [10] and therefore useful only at large internuclear 
distances. Another complication is that for the homonuclear case, every gerade energy E g (wave function 
symmetric under exchange of nuclei) has a counterpart ungerade solution (wave function antisymmetric 
under exchange of nuclei) whose energy E u has exactly the same 1/R expansion. This makes the calcu- 
lation of exchange energy splittings AE = E u — E g very elusive to calculate at large R, although there 
are specialized methods for recovering these splittings (e.g. see [11]). 

Even recently, there has been examination of series in small R limited to the ground state short-range in- 
teraction energy [12] but we still have no further insight into the actual mathematical nature governing the 
energy eigenvalues. We also cite the work of Demkov et al. [13] but their analytical solutions correspond 
to a peculiar charge ratio depending on the internuclear distance and therefore not physically useful. 

Thus, complete analytical solutions of the eigenstates of H^, in areas of molecular interest, such as e.g. 
the region near the equilibrium internuclear distance (bond length) of the ground state remain elusive. 

A wide variety of numerical methods have been used to solve the H^" problem in this case. For example, 
Bates, Ledsham and Stewart [5] used recursion and continued fractions. Hunter and Pritchard [6] used 
matrix methods and Rayleigh quotient iteration. Madsen and Peek [7] used power series and associated 
Legendre expansions to set up two equations whose simultaneous solution then gave the two eigenpa- 
rameters. An accurate way to obtain energies and wavefunctions for the one-electron two-center problem 
is provided by the program ODKIL conceived by Aubert-Frecon et al. [14, 15] based on a method by 
Killingbeck. As of the 1980s, it was possible to calculate the eigenenergies and the eigenf unctions of 
the discrete states of H^" with a rapid FORTRAN program. Yet, complete analytical solutions have so 
far remained elusive: the classical iV-body problem cannot be solved in closed form for N > 3 and the 
quantum counterpart is even worse by virtue of being an eigenvalue problem. 

The approach used here is called "experimental mathematics", an unorthodox approach involving multi- 
disciplinary activities by which to find new mathematical patterns and conjectures. The goal in this context 
is to search and find mathematical structures and patterns to be re-examined with more "rigor" at a later 
stage. The level of rigor is of course relative: in dealing with a difficult problem in applied Mathematics, 
we cannot approach the level of rigor demanded in number theory. Nonetheless, we desire demonstrations 
sufficiently convincing to the molecular physicist. 

The present work will involve a combination of methods, results and procedures from different areas. 
We first start with results from what is called: dimensional scaling. It has been known for some time 
that the Schrodinger wave equation can be generalized to an arbitrary number of dimensions D which 
can be subsequently treated as continuous variable [17, 18]. In the limit as D — > 1 + , the hydrogen 
molecular ion becomes the double well Dirac Delta function model which can be solved exactly [25] in 
terms of the Lambert W function [19,20]. Dimensional scaling applied to H^" has been studied at length 
by Hershbach's group [21-24], in particular, by Frantz [21], Loeser and Lopez-Cabrera [22,23]. The latter 
work provides even more insight into the mathematical relationship between the real at D = 3 and its 
one-dimensional limit. 

Next, armed with the information provided by dimensional scaling, we will return to the real three- 
dimensional formulation of Aubert et al. [15]. This formulation is re-examined using a Computer Algebra 
System (CAS) within the approach of experimental mathematics: patterns and results are obtained. The 
CAS used is Maple because it is readily available to us but the results could also be implemented on other 
systems. The resulting series expansions are verified numerically and algebraically. In particular, we will 
demonstrate that our results are independent of choice of basis and basis size and consequently completely 
general. The end-result will be then analytically compared with the one-dimensional result and put on a 



near equal footing allowing us to find the mathematical category to which belong the eigenvalues of . 
In view of the type of solution obtained, a tentative "physical" picture is associated with the analytical 
solutions. A summary with concluding remarks is made at the end. 



2 Preliminaries - Dimensional Scaling 

The D — ► 1 + version of H^" [17, 18] is given by the double Dirac delta function model: 

where Za = q and Zb = X q. The ansatz for the solution has been known since the work of Frost [26]: 

$ = Ae~ d \ x \ + Be~ d \ x - R \ (2) 
Matching of ip at the peaks of the Dirac delta functions positioned at x = 0, R when (A = 1) yields: 



q — d qe dR 
qe~ dR q - d 



(3) 



and the energies are thus given by: 

E± = -d 2 ±/2 where d± = q[l ± e~ d±R ] (4) 

Although, the above has been known for more than half a century, it was not until Scott et al. [25] that the 
solution for d± was exactly found to be: 

d± = q + W(±qRe~ qR ) /R (5) 

where ± represent respectively the symmetric or gerade solution and the anti-symmetric or ungerade solu- 
tion and W is the Lambert W function satisfying W(t)e w ^ = t [19,20]. This function first introduced by 
Johann Heinrich Lambert (1728-1777), a contemporary of Euler, has been "invented" and "re-invented" 
at various periods in history but its ubiquitous nature was not fully realized within the last decade or so. 

For example, the W function appears in Wien's Displacement Law of Blackbody radiation. In general, it 
has appeared in electrostatics, statistical mechanics, general relativity, radiative transfer, quantum chromo- 
dynamics, combinatorial number theory, fuel consumption and population growth etc. . . (e.g. see ref. [27] 
and references herein). 

More recently, the Lambert W function has also appeared in "linear" gravity two-body problem [28] as 
a solution to the Einstein Field equations with one spatial dimension and one time dimension (1 + 1). 
The present work also includes a generalization of the W function. Recent work [30] shows that the W 
function can be further generalized to express solutions to transcendental algebraic equations of the form: 

P N (x) 

exp(±ca;) = (6) 
Qm{x) 

where Pn{x) and Qm(x) are polynomials in x of respectively degrees N and M and c is a constant. The 
standard W function applies for cases when N = 1 and M = and expresses solutions for the case of 
equal charges for eq. Q or equivalently the case of equal masses for the two-body 1 + 1 linear gravity 
problem. The case of unequal charges or unequal masses corresponds to cases of higher N and M values. 
This form also expresses a subset of the solutions to the three-bo&y linear gravity problem [29, 30] where 
one deals with transcendental equations of the form (0 where M, N — > oo. 



Some insight into the mathematical nature of the eigenenergies of H^" is revealed by the fact that the 
eigensolutions for the electronic energies at D — ► 1 + and D — > oo actually bound the D = 3 ground 
state eigenenergy of [21, 22] as shown in Figure 1. Moreover, the latter can be estimated by a linear 
interpolation formula [23] : 

£ 3 (i?)~^i(f) + |s„(^) (7) 

This formula agrees with the numerically accurate eigenenergy (as given by program ODKIL or the work 
of D. Frantz) to within about 2 or 3 digits for the range of R near the bond length. The result at D — > oo 
involves the extrema of a Hamiltonian expression [23, eq.(58)]. We re-examined this result. One has to 
consider a region of R divided by R c = |v3- For R < R c , the root is determined by the root of a 
quartic polynomial [23] and the result for R > R c is determined by a sixth degree polynomial. Thus, the 
result at D — ► oo is algebraic. On the other hand, the result at D —> 1 + is in terms of an implicit special 
function, which is the Lambert W function. Given how well this interpolation formulation works, this 
already suggests what is the mathematical nature of the eigenenergies of the true hydrogen molecular ion 
(D = 3). 

We can state this in view of the work of Frantz et al. [21] who showed that the D-dimensional problem 
could be decoupled into two coupled ODEs for 2 < D < oo and how a particular energy eigenvalue for a 
given D could exactly express the solution of another eigenvalue for an excited state at a dimension D + 2 
through a precise re-scaling. 



3 Three-dimensional 



3.1 Starting Formulation 



The Schrodinger Wave Equation for in atomic units is given by: 



ijj = 



(rA + r B )/R , 1 < e < oo 
(r A -r B )/R , -1 < t; < 1 



. r a r B 

As mentioned before, this is separable into prolate-spheroidal coordinates: 

V 

< (f> < 2vr 

Qi = R(Z A -Z B ) 
Q 2 = R{Z A + Z B ) 

We can write the ansatz for the eigensolution: 

^, V A) = MO M (r], 4>) = A(C)G(r ? )e ± ^ 
which allow us to obtain two coupled ODEs: 



— (a- 2 )— 

drj \ ^ drj 



-((f-l)- 



m 



1-7? 2 

9 

m 



+ p 2 r] 2 -QiV-A 
v 2 t 2 + Q 2 t + A 



M(r],(j)) = 
A(£) = 



where A is the separation constant and the eigenenergy E is expressed as: 

R 2 



(8) 



(9) 



(10) 



(11) 



Note that: 



lim A 

lim p 2 

R^O 







(12) 
(13) 



Although the set of quantum numbers (n, i, m) - the united atom quantum numbers - can be used to 
identify the eigenstates, as is the case for e.g. program ODKIL, it must be emphasized that only the 
magnetic quantum number m is a good quantum number (resulting from the azimuthal symmetry of 
about its internuclear axis). 

We follow the treatment of Aubert et al. [15] and consider the following basis expansion for the rj coordi- 
nate: 

M{r 1 A)= ]T/rn m (^) (14) 

k=m 

where Y" fc m are the Spherical Harmonics. Injection of the above basis into the ODE governing M in rj 
leads to the creation of a symmetric matrix T whose determinant must vanish when p and A satisfy the 
eigenvalue problem: 

[F(p,A)]\f\=0 (15) 



where 



Fi,i(p,A) 
F i+ i,i(p, A) 
F i+2>i (p, A) 



+ 1) + p 
-RQi 



i 2i 2 - 2m 2 + 2i - 1 



(2i + 3)(2»- 1) 
i + m + l)(z — m + l" 1 ^ 



A, 



(2i + 3) 



(2i + l)(2i + 3) 
i + m + l)(i - m + + m + 2)(i - m + 2)^ 1/2 



(2z + l)(2z + 5) 



are the non-vanishing matrix elements of T . If Q\ = Zj± — Z% = i.e. the homonuclear case, then the 
pentadiagonal matrix T divides in even and odd tridiagonal matrices in terms of x where x = p 2 with no 
explicit dependence on the internuclear distance R (although this is not true for the other ODE in £). For 
the £ coordinate, we use a basis of Hylleraas functions, i.e. in terms of Laguerre polynomials: 



A(o = e -f«- i )[2p(e-i)r/2 y, q 

n=m/2 



pm, 

n-(m/2) z -n-(m/2) 



[Mt - 1)] 



(16) 



where 



[y(p,A)]\c\ = o 

[y(p,A)} = [Q(p,A)] +pm 2 [B(p)r l 



(17) 



and 



h,i(p) = 4p + 2i + l, 



772 777. 

bi,i+i(p) = b i+1>i (p) = - + + _ + 



r iji+ i(p, A) 



C* + l)|'^-« 



1/2 



| + ^- + i + i? Q 2 - P 2 + A, 



1/2 



RQ2 . 

X Z 

2p 



When m = 0, the matrix is tridiagonal. For m ^ 0, one has to consider the inverse of the matrix £>, which 
is not a band matrix. 

Of course, we realize that this choice of basis is only one of several possible choices. The results obtained 
are valid provided the results are independent of the size of the basis and the choice of basis. 



3.2 Recurrence Relations 



The following relations apply to the homonuclear case and when m = in which case, the band matrices 
are purely tridiagonal matrices. These are governed by recurrence relations namely (Al) and (A.2) of 
reference [15]: 



detpfto] 
det[9Jti] 
det[Wl k ] 



1 

mil 



m k ,k detpJVi] - m k -i,k m fcifc _i det[9Jt fc _ 2 ] 
Thus for m = 0, we have the following: 

det(Jo) 
detCVi) 



1 

R 

P 



l-2p + 2R-p z + A 



R 



(18) 



;2fc + l) I — - & - 1 -2p J + k + 2R-p z + A) det(y k ) 



For the even I case, we have: 
det(^ r eo) 
det(.Fei) 

det(J"e fc+ i) 



-A 

. p 2 (8 k 2 + 4 k - l) \ 
-2M2 t+ l) + ^_ 1)(4t + 3 ; -Ajdet(^) 

. p A (2k- l) 2 A; 2 det(^e fe _i) 



(4k- l) 2 (4k-3) (4fc + l) 
Denning det(JP"o,j) = det(jF(i)) for the odd I case, we have: 



(19) 



det(J r oo) 
det^oi) 

det(J"o fe+ i) 



1 



5 



-(2fc + l)(2A; + 2) + 



p 2 (2 (2/c + l) 2 + 1 + 4k 
(4k + 5) (4k + 1) 



A det(JT 0A .) 



p 4 fc i (2fc + l) / det^Qfc,!) 
(4fc + l) 2 (4k - 1) (4fc + 3) 



(20) 



Note that the radial equations for the hydrogen atom are governed by two-term recurrence relations. Thus, 
it suffices to find an eigenenergy such that the coefficient a k+ \ of the basis of Laguerre functions is zero. 
This in effect truncates the infinite series into a polynomial and consequently closed form solutions for 



the eigenstates are obtained of the hydrogen atom. This is not possible for which is governed by 
three-term recurrence relations no matter what the choice of basis. 



The band matrices for and their determinants have been injected into a computer algebra system. The 
determinants det(J^) and det(J r j) (even or odd) for i = 1, 2, 3 . . . are multivariate polynomial-like in A 
and p. The determinants det(^ r j) are true polynomials in A and p 2 . On the other hand, although det (3^) is 
a polynomial in A, it has also negative powers for p and thus akin to a Laurent series (Laurent polynomial) 
in p. 

It is possible to eliminate one of the unknowns by obtaining a resultant of the two determinants det (3^) 
and det(J r j). If a and b are polynomials over an integral domain, where et 2 (rational) polynomial equa- 
tions in 2 unknowns A and p. 

n m 

a = a n Y[(x - an) b = b m JJ(a; - /%) 

i=l i=l 

Then 

n m 

resultant(a, b, x) = a^b^ J — Pj) 

i=l j=i 

This can be computed from the Euclidean algorithm or determinant of a Sylvester matrix and its roots will 
be common to those satisfying the original set of polynomials. Since both expressions are true polynomials 
in A only, the resultant must be in A. E.g. for i = 2 (i.e. 2x2 matrices) 

resultant (det (3^ ), det(.Fej), A) = 

64/1225 p 8 + 512/245 p 7 - 256/245 (R - 27) p 6 - 128/245 (56 R - 369) p 5 
+ 64/245 (2911 - 1037 R + 27 R 2 ) p 4 + 32/245 (13580 - 9405 R + 948 R 2 ) p 3 
- 32/245 (140 R 3 - 17780 + 24010 R - 5571 R 2 ) p 2 - 64/7 (20 R 3 - 224 R 2 
+ 481 i?- 161) p- 4128/7 R 3 + 19968/7 R 2 + 16 i? 4 - 2880 R + 304 
+ 16/7 R (28 ,R 3 - 347 R 2 + 791 R - 252) /p + 4 R 2 (20 R 2 -W8R + 85) /p 2 
+ 8 R 3 (4 R - 9)/p 3 + 4 R A jp A 

i.e. a Laurent polynomial in p with coefficients in R only. When the size of the i x i increases, the size 
of the resulting expression increases dramatically (expression swell). However, from a numerical point of 
view, the most useful outcome comes from numerically solving the simultaneous expressions for det (3^) 
and det(J r j) since i must be sufficiently large to give a sufficiently good result near the bond length. In 
Maple, this can be done using the f solve procedure. To find the minimum energy for the ground state, 
it is a matter of getting derivatives of these determinants with respect to R. Combining the latter with the 
condition: 

where E T = E elec + 1/R (21) 



OR 

we get five equations in the five unknowns R, A, p, ^ and J^r. The result has been calculated using a 



small Maple program. In atomic units, these are: 

R = 1.997193319969992... 
Eminimum = -0.6026346191065398... 

Note that the electronic energy, evaluated at R = 2.0 a.u. for comparison, is as expected exactly the 
reference tabulated value of Madsen and Peek [7] i.e. —0.6026342144949. An indirect way of ascertain- 
ing the accuracy of electronic energies is to use these values in an adiabatic standard scheme to obtain 
vibrational energies which are directly comparable to highly accurate values provided by approaches that 
do not involve the separability of the electronic and nuclear motions (e.g. [31-33] and [34,35]). This has 
been done [36] and comparisons with values from the literature are displayed in table 13 .21 



Table 1 : Ground State Vibrational Energies 

Energies are in a.u., differences in cm -1 (1 a.u. = 219474.63 cm -1 ) 
a) ref. [31], b) ref. [32], c) ref. [33], d) ref. [34], e) ref. [35] 





Quantum 


Present 




Differences 


System 


Vibrational 


Adiabatic 


Literature 


AE 




Number v 


Values 


Values 


(cm" 1 ) 


H+ 





-0.597138471 


-0.597139055 a) 


-0.13 








-0.5971 39063 123 6 ) 


-0.13 




1 


-0.587154167 


-0.587155679212 6 ) 


-0.33 


D+ 





-0.598788594 


-0.5987876(1 l) a ) 


+0.22 








-0.59878878433 l c ) 


-0.04 


HD+ 





-0.597897521 


-0.5978979685 d ) 


-0.10 








-0.5978979686 e > 


-0.10 



In fact, given how heavy the nuclear centers are with respect to the electron, clamping the nuclear centers 
is a very good approximation for the quantum three-body problem represented by with the following 
caveat: the approximation that the nuclei are clamped fixed in space creates a symmetry under exchange 
of nuclei in the homonuclear case. A different picture arises when the movement of nuclei is considered. 
The mere movement of the nuclei breaks the symmetry under exchange of nuclei and thereby leads to a 
localization of the states. In this case, the work of Esry and Sadeghpour is instructive [37]. 

However, if one stopped here, there is no pattern from an analytical point of view. E.g. setting x = p 2 
and examining det(^ r ej) at low order in A, we have: 



at i = 2: 



at i = 3: 



at i = 4: 



l/35p 2 (-70 + 3p 2 ) + (6 - 6/7p 2 )A + O (A 2 ) 



JL p 2 (1848 - 126 p 2 + /) + (-120 + ^ p 2 - A p 4 )A + O (A 2 ) 



p 2 (-2162160 + 173316p 2 - 2772 p 4 + 7p G ) 



1287 

+(5040 - 1032/ + ^P A ~ ^P 6 )A + O (A 2 ) 



If we look at A = and grab the leading coefficient p 2 , we have the sequence 
—70, 1848, —2162160 .... Not only are the coefficients increasing dramatically in size, they also alter- 
nate in sign. Although the roots A,p of these determinants det(3^i) and det(J r j) converge with increasing 
i, the actual coefficients of these determinants and especially those of the resultant increase in size becom- 
ing more and more cumbersome although a CAS can handle them (up to a point). 



Moreover, we have made a particular choice of basis and the combined set of polynomial-like expressions 
for the determinants though numerically useful could be viewed more as a computational "model" rather 
than anything truly representative of the wave function. If we stop here, we see no pattern. Insight comes 
from inverting the problem. 



3.3 Roots of Determinants 

The three-term recurrence relations for det(3^) or det(^ r j) cannot be solved in closed form. We start with 
det(J r j) because it is easier and has no explicit dependence on R. Upon careful scrutiny of eqs. (f!9l and 
(120b . the term in det^k-i) has a coefficient in p 4 whereas the term in det^k) has terms at order p 2 . Let 
us assume that p is small, which is indeed the case for small R. We can therefore neglect the last term in 
det(J r / c _i) and the resulting two-term recurrence relation becomes trivial to solve. It is merely a matter 
of compounding the multiplicative terms of the recursion: 

det(^) « ( -i/n + 1) + a - < 22 > 

Solving for A such that det(.Fefc) = yields: 

a = -a w-wir + y +o{p,) <23> 

We can clearly identify the R — > limit with I = 2 j. Similarly, for the odd case, we have: 

» (-i)'n ( 2 o+^+i) + A - %tw'j% ) <24) 

A=- W + l )W+ Q + $±^f + 0<S) (25, 

We can clearly identify the R — > limit with I = 2 j + 1. Thus, although £ is only a valid quantum number 
in the united atom limit, it is nonetheless feasible to use it to identify an eigenstate as an expansion for 
small p (and small R). 

By the implicit function theorem, det(.Fj) = =4> A = A(p 2 ). Moreover, the structure of the recurrence 
relations for det(J r j) and det(3 ; i) namely eqs. ( fT9l . d20t and (Tf8b tell us that all these quantities are i th 
degree polynomials in A. If one can find all the values of A such that these determinants are zero, the 
latter are clearly known by the fundamental theorem of algebra. If det(J r j) as a formal series in x where 
x = p 2 , we can use reversion of power series to obtain an analytical solution. This is the best possible 
analytical result. E.g. , we consider: 

x = cos(x) =>■ x/cos(x) = A where A = 1 (26) 



^ = A--A' + -A°+ ... 

The reverted series of x in terms of A can be obtained in a number of ways including Lagrange's method [4] 
and represents the best possible representation of an analytical solution to the root of eq. d2^l . Formally, 
the infinite series in A is a complete solution. The issue of getting numbers for e.g. A = 1 is a matter of a 



summation technique. Solutions by reversion of power series are possible via Maple's solve command. 
E.g. inverting det(.Fe3) yields: 

l/3x + —a? + 0(a?), -6 + — x - — x 2 + O (x 3 ) , 
1 135 K ' 21 9261 v ' 

39 8 o „ / q\ 

-20 + -X-— x 2 + 0(x 3 ) 

where x = p 2 . To first order in x (or p 2 ), we recover the solutions in eq. d23l for respectively £ = 0, 2, 4. 
The action of inverting det(x~ej) produces i solutions to order 0{x % ). E.g. if we isolate the I = solution 
obtained from det(x~ej) for i = 3, 4, 5, 6, we obtain: 

i = 3: l/3x + -|-x 2 + 0(x 3 ) 
135 

j = 4: l/!«4« ! ^^0(«') 

i = 5: l/3x + x 2 + x 3 x 4 + O (x 5 ) (27) 

7 135 8505 1913625 v ; 

2 2 4 o 26 4 92 = n , BN 

, _ 6 : 1/3. + x» + *» - -^tOM 



2 2 4 o 26 4 92 = 

i — > oo : l/3xH x H x x x +... 

7 135 8505 1913625 37889775 

What is important to note is that the coefficients are stablel Letting i —> i + 1 adds a term of order 0{x t+l ) 
to the series and yields an extra solution for £ = 2 (i + 1). By re-injection of this solution to within order 
0{x % ) into det(x~ej) with computer algebra, one can see that det(x~ej) is satisfied, term by term to within 
that same order. Conversely, the coefficients of A for a particular choice of £ even can be obtained from 
this simple algorithm: 

1 . Select value of j and £ = 2 j and desired order N. 

2. Set ao and a\ according to eq. <I23I > 

a = -£{£ + 1) 
ai : 



8j 2 + 4j - 1 



(4j-l)(4j + 3) 

3. For i = 1 to (N - 1) 

(a) Let A tria i = Ylk=o a k xk + a «+i Note that aj+i is symbolic and not yet determined. 

(b) Substitute A tr i a i into det(x _ e^ + j + i). 

(c) Isolate coefficient for x t+l . 

(d) Solve for a^+i such that this coefficient is zero. 

A counterpart result also holds for the odd case of £ i.e. for det(^ej). This simple algorithm allows us to 
yield the series solution for A for any given choice of £. At the same time, the solution of this algorithm 
implies that det(x~ei) = is formally solved. 

It must be emphasized that increasing i merely means adding basis functions. There are no singularities 
between the two nuclei of H^", and we can expect the wave function to be not only continuous but also 



continuously differentiable in that regime i.e. we expect no surprises with the basis functions as i — > oo. 
As the estimates for A and p are closer and closer to the true values of the eigenparameters, the magnitude 
of the coefficients fi of eq. (fT3l become smaller and smaller as i — > oo. In this limit, the basis set is a 
valid representation of the true wave function. 

The first 10 coefficients of the series for A(x) where x = p 2 for £ = are: 

2 n 4 o 26 4 92 . 

A(x) = 1/3 x + x 2 + x 3 x 4 x 5 (28) 

V ' ' 135 8505 1913625 37889775 V 

513988 fi 122264 7 57430742 8 
■ x° + —— x' + — — — , „„, x s 



9050920003125 11636897146875 62315584221515625 

26237052532 ^ 1550889714543116 10 



1566426840576238265625 213229853673440433908203125 

and our computer algebra programs allow us to generate many more such coefficients. The first three non- 
vanishing terms of the Taylor series for A(p 2 ) have already been published for cases of small p consistent 
with small internuclear distances R [38-40]. We now claim that the present algorithm provides a means 
of generating the Taylor series of A in small x where x = p 2 , the result being valid as i — > oo and thus 
independent of the size of the truncated basis. Later, we will demonstrate it to be independent of the actual 
choice of basis. However, the first test concerns numerical vindication. 



3.4 Numerical vindication of the Series for A(p 2 ) 

To vindicate the series, we obtain data entries of R, p and A from program ODKIL and inject the data 
entries of p into the series solutions for A. We then compare the latter with the value of A obtained 
from ODKIL for a given state. This is done for the ground state and a few excited states as shown in the 
following tables. The results for the ground state i.e. Is <r g (n = 1, 1 = 0, m = 0) are shown in tableland 
those of state 2s a g (n = 2, I = 0, m = 0) are shown in table fj] demonstrating that the same series of A 
for a given £ works for more than one state. The results for the excited state 2p a u {n = 2, £ = 1, m = 0) 
vindicate the series solution for t = 1. The results for state 3da g (n = 3,1 = 2, m = 0) vindicate the 
series solution for A for £ = 2. 

In all cases, we can see that the series obtained for A{p 2 ) works indeed like a Taylor series, working 
very well for small p. Beyond a certain value of R, the series solution rapidly degenerates. Nonetheless, 
e.g. for the ground state, the series solution works well near the bond length (around R = 2 which is 
underlined) and beyond. Degradation of the series becomes apparent at R = 5. 

The question arises as to whether or not the series coefficients of A{p 2 ) follow a pattern. We have found 
none so far. The pattern of the changing signs +, — is not one of alternating series and thus this function 
is unlike all the special functions known in the literature (such as e.g. [41]). 

Nonetheless, there is something of a pattern for a given series when modifying the quantum number £, 
term by term. The first two terms ao and a\ follow a pattern in I according to e.g. d23l for even £. No 
such simple pattern exists for the next term a^,. However, if one solves for 02 in terms of ao and a\ for 
a high value of £, say £ = £ ma x, one obtains a polynomial formula for 02- If one then substitutes the 
general formulae in £ for ao and a\ into this polynomial expression for a^. it will correctly generate the 
coefficients 02 not only for £ ma x but for all £ = 0, 1, 2 . . . £ m ax- At some point, the resulting formula 
will break down for a value of £ > £ ma x- This "triangular" relationship, - useful because one often does 
calculations within for a limited range of £ - indicates that: 



Table 2: Convergence of "Taylor series" of A(p 2 ) 
Ground State: Is a g (n = l,£ = 0,m = 0) 



1 F? 

1 .ft 


ODKIL (accurate) 


series yzxj lerms ) 


P 


A 


A 














0.5 


.46569679 


.729927345e-l 


.7299273577e-l 


1 A 

i.U 


.851993637 


.249946241 


.249946 15 /4 


1.5 


1.18537488 


.498858904 


.4988725127 


]Z0 


1.48501462 


.811729585 


.8118596153 


2.5 


1.7622992 


1.19023518 


1.190951531 


3.0 


2.02460685 


1.64100244 


1.643819599 


4.0 


2.52362419 


2.79958876 


2.822919217 


5.0 


3.00919486 


4.37769375 


4.491055954 


10.0 


5.47986646 


20.1332932 


25.05609231 


20.0 


10.4882244 


90.0528912 


-1147.477000 



Table 3: Convergence of "Taylor series" of A(p 2 ) 
State: 2s a ( n = 2,£ = 0,m = 0) 



R 


ODKIL (accurate) 


series (20 terms) 


P 


A 


A 














0.5 


0.241110452 


0.194282436e-l 


0. 194282436 le-1 


1.0 


0.459850296 


0.71 1543 142e-l 


0.7115431427e-l 


2.0 


0.849546791 


0.248466171 


0.2484661714 


3.0 


1.19791141 


0.510154273 


0.5101542740 


4.0 


1.51924947 


0.8535318 


0.8535318053 


5.0 


1.82176362 


1.28400188 


1.284001886 


10 


3.19930169 


5.12935962 


5.127249696 


15 


4.51129751 


12.4337232 


-17315.20146 



which places us beyond eq. d23t (or d25l ) which determine ao and a\ only. However, this is subject of 
further exploration elsewhere. 

The range of the series solution can be considerably improved by modifying the recurrence relation for 
det(J r e) like so: 



det(J r eo) 
det(.Fei) 

det(J"e fc+ i) 



1 

y 

3 



y --A 



y(8k 2 -l + 4k) \ 
-2k(2k + l) + ^ 1 - A] det(Fe k ) 



(4k- 1) (4k + 3) 

xy(2k- l) 2 k 2 
(4k - l) 2 (4k - 3) (4k + l) 



(29) 
(30) 

(31) 



det^efc-i) 



where it is understood x = y = p 2 but it is only x which is treated as a perturbation. This is simply a 
different representation denoted A = A(x,y) but which represents the same function A(p 2 ). Modifying 



Table 4: Convergence of "Taylor series" of A(p 2 ) 
State: 2pa u (n = 2,t = 1, m = 0) 



R 


ODKIL (accurate) 


series (70 tprmO 


p 


A 


A 


U 


u 


-z 


-L 




.254186316 


-1.96120498 


-1 961204981 


1.0 


.53141962 


-1.83001042 


-1.830010419 


2.0 


1.15545177 


-1.18688939 


-1.186889387 


4.0 


2.35889913 


1.53846448 


1.538464473 


6.0 


3.43970785 


5.92793017 


5.927930398 


8.0 


4.4671459 


12.0646853 


12.07439611 


9.0 


4.97308004 


15.8356448 


16.60977070 


10 


5.476774 


20.0920989 


58.89905749 


20 


10.4882239 


90.0528776 


.7649129703el3 



Table 5: Convergence of "Taylor series" of A(p 2 ) 
State: 3d cr g (n = 3, £ = 2, m = 0) 



R 


ODKIL (accurate) 


series (20 terms) 


P 


A 


A 




0.5 
1.0 
2.0 
4.0 
6.0 
8.0 
10 
20 


0.166934253 
0.335547827 
0.686698811 
1.51188304 
2.37168861 
3.09069127 
3.69538523 
6.12806789 


-6 

-5.98541087 
-5.94115241 
-5.75530105 
-4.86085811 
-3.43229937 
-2.07684281 
-0.874720469 
7.31365225 


-6. 

-5.985410869 
-5.941152409 
-5.755301048 
-4.860858108 
-3.432299419 
-2.076688125 
2.071237971 
5232651466. 



slightly our previous algorithm, we obtain e.g. a modified series solution for A(x, y) for t = 0: 

A = y.-— yX + — y2 ^ (32) 

3 15 (2 y- 63) 375 (2 y - 63) 3 (2 y - 231) 

28 y 3 P 3 (y) x 3 



121875 (2y-63) 5 (2 y - 231) 2 (2 y - 495) 
where the polynomials Pk(y) of order k are given by: 

P^ y ) = 94 y- 44121 

P 3 ( y ) = 166376 y 3 + 16398492 y 2 + 131745081006 y - 13685763372435 

Note that if we inject y = x into the above and make a Taylor series expansion in x, we simply recover 
the series solution in x = p 2 obtained in eq. (|2"71 for £ = 0. Since the radius of convergence is determined 
by the closest singularity or branch point in the complex plane, we have 



2 y - 63 = (2 p 2 - 63) = p « 5.6 



We note that the sequence of numbers 63, 231, 495, 855,. . . which appear in the denominator have a pattern 
which can be found using the gf un package [45]. This demonstrates that these numbers fit a holonomic 
function and it is found that these fit the pattern: 

3(4j + 3)(4j + 7) (33) 

We recognize it as one of the terms which appear in the recursion relations for det(.Fefc) 
i.e. (4k — l)(4fc + 3) with k = j + 1. However, no pattern has (so far) been found for the polyno- 
mials Pk{y)- Nonetheless, our computer algebra routines allow us to generate this series to relatively high 
order. 

Next, the sum can be calculated using non-linear transformations known as the Levin or Sidi transfor- 
mations. The latter involves a series transformation by which one can accelerate the convergence of 
a series and even sum divergent series (e.g. see the work of [43, 44]). We take the point of view 
that a Taylor or asymptotic series has all the desired "information", getting numbers from the series 
is a matter of a summation technique. These transformations are available in the Maple system as 
Nonlinear Trans format ions. 

The best results for the ground state are obtained by applying a Sidi d transformation in x compounded 
with y as shown in table |6] Even when the modified series behaves badly, the result from the Sidi d 
transformation provides reliable numbers. The results hold up remarkably well all the way up to R = 10 
and beyond. Beyond R = 10, the asymptotic series expansions as e.g. listed by Cfzek et al. [9] are more 
useful. What is important in our case, is that our series solution works so well around the bond length and 
the intermediate regime. 



Table 6: Convergence of Series A(x,y) 

Ground State Revisited: Is cr (n = 1,£ = 0, m = 0) 



R 


A 


series (12 terms) 


ODKIL (accurate) 


Sidi-d 


1 


0.24994624090 


0.2499462409 


0.2499462410 


2 


0.81172958404 


0.8117295840 


0.8117295850 


3 


1.6410024366 


1.6410024369 


1.6410024370 


4 


2.7995666114 


2.7995887586 


2.7995887590 


5 


3.9638237398 


4.3776938960 


4.3776937530 


6 


-4.6313683166e+03 


6.4536051398 


6.4536037430 


8 


-3.7137673759e+12 


12.2262006172 


12.2261746150 


10 


-1.5608159299e+33 


20.1339450995 


20.1332931780 


15 


1.005460041 le+15 


48.8656127918 


48.8223535290 



3.5 Solution for A(R, p) 

Although we have eliminated one of the unknowns i.e. found A(p 2 ) such that the determinantal conditions 
for det(J r j) are satisfied, there is still the remaining determinant det(3^) to address. The recurrence 
relations for det(3^) of eq. (Tf8l) depend on the internuclear distance R and have more structure than those 
of det(^ r ej) of eq. ( fT9b or det(J r Oj) of eq. d20t . Nonetheless, we proceed in parallel to what we did for 

A( P 2 ). 

To start with, we ignore the term det(3 ; fc-i) and solve the resulting two-term recurrence relation since all 
linear recurrence relations of this type are solvable in terms of the roots of the characteristic polynomial 



obtained by assuming a det(3 ; i) = f l and then solving for /: 



Hptn , x _ ( _ ,k r((2 k P + y + + x/2 P )) r((2 fe p + y + - x)/(2p)) 

Wj " ( j r(-(y_-jr)/(2p))r(-(r- + x)/(2p)) (34) 



where 



X = v / 2p i -p 2 + R 2 + 2 Ap 2 
Y + = 2p 2 + p - R 

y_ = 2 P 2 - P + r 

and T is the Gamma function [41]. This result bears some resemblance with the outcome of solving the 
eigenvalue problem for the hydrogen atom. In this case, solutions to the ODE for the radial equation in 
the radius r can be expressed in terms of hypergeometric functions. Matching the asymptotic solution at 
r — > oo with the regular solution at r — > necessitates the elimination of the irregular solution by forcing 
one of its coefficients - also expressed in terms of the Gamma function - to be zero (e.g. see [42]). In our 
case (as in the case of the hydrogen atom), it is a matter of ensuring that the arguments for one (or both) 
of the Gamma functions in the denominator of the expression above to be — j where j = 0, 1, 2 . . . Thus, 
solving for A, we find that: 

A(R,p) « p 2 + 2(l + 2j)p + 1-2R + 2j + 2j 2 - ^ii±M . ( 35 ) 

P 

What remains is the identification of j. Next we treat term det(3^_i) as a perturbation formally by 
multiplying it by A with the understanding that (A = 1). For j = 0, the series solution for A(R,p) is: 

A(R,p) = (p+l) 2 -2R-* - - ^V/^T m 06) 

p 2 p (2p z + 2 p — R) 

(p-R) 2 P i (R,p)X 2 



(8 p (2 p + 2 p 2 - R) 3 (3 p + 2 p 2 - R)) 

(p-R) 2 P w (R,p)\ 3 

(16 p (2p + 2p 2 - R) 5 (3p + 2p 2 - R) 2 (4p + 2p 2 - R)) 

where 

P 4 (R,p) = Up 4 + (13- 12R) p 3 + R(2R- 17) p 2 + 7 pR 2 - R 3 
P w (R,p) = 584p 10 - 8(116 R- 233) p 9 + 2(256 R 2 + 969 - 1878 R)p 8 

- 4 (28 R 3 - 722 R 2 + 1184 R - 165) p 7 
+ R (8 R 3 + 4678 R - 1056 R 2 - 1897) p 6 
+2R 2 (92 R 2 - 1200 R + 1145) p 5 

- R 3 (12 i? 2 - 678 R + 1513) p 4 - 2 i^(50 R - 297) p 3 
+ R 5 (6 J R-139)p 2 + 18i? 6 p - R 7 

The series looks complicated and the presence of singularities at every —R+i p+2 p 2 = for i = 2, 3, . . . 
already tell us that this function is unlike most special functions in the literature. However, the series gives 
very good results as shown in table with only 4 terms. It does not need any convergence acceleration 
summation methods at large R. The results of Table[8]for state 2s a g (n = 2, i = 0, m = 0) show us that 
A(R,p) works well for large R but diverges for small R. Also shown in the table are the results of the 
Sidi d transformation which considerably improves the series solution for small R. 

What remains is to identify the meaning of the number j. By checking the solution for excited states, we 
find out empirically that: 

j = n - £ - 1 (37) 



where n is the united atom quantum number. This number j is a valid quantum number for the separated 
atom limit [8, eq.24,p.666]. Thus, just as we match the outward and inward radial solutions for the radial 
ODE for the hydrogen atom by which to determine the eigenvalue, the eigensolution for results from 
matching A(p 2 ) governed by the united atom quantum number £ with A(R,p) governed by the separated 
atom quantum number j = n — £ — 1. 



Table 7: Convergence of A(R, p) 

Ground State: lsa g (n = 1, 1 = 0, m = 0) 



R 


ODKIL (accurate) 


series (4 terms) 


P 


A 


A 


0.5 


0.46569679 


0.729927345e-l 


0.7299778055e-l 


1.0 


0.851993637 


0.249946241 


0.2499480309 


ZO 


1.48501462 


0.811729585 


0.8117297560 


5.0 


3.00919486 


4.37769375 


4.377693772 


10 


5.47986646 


20.1332932 


20.13329314 


20 


10.4882244 


90.0528912 


90.05289034 


30 


15.4919739 


210.034597 


210.0345960 


40 


20.4939187 


380.025707 


380.0257060 


50 


25.4951064 


600.020452 


600.0204512 



Table 8: Convergence of A(R, p)State: 2s a g (?i = 2, £ = 0,m = 0) 



R 


ODKIL (accurate) 


series (5 terms) 


Sidi-d 


P 


A 


A 


A 















0.5 


0.241110452 


0.0194282436 


-7.5860659805e+02 


0.0211395500 


1.0 


0.459850296 


0.0711543142 


-2.1405738045 


0.0718499907 


2.0 


0.849546791 


0.2484661710 


0.23871213039 


0.2485999772 


3.0 


1.197911410 


0.5101542730 


0.50971077859 


0.5101743643 


4.0 


1.519249470 


0.853531800 


0.85348392428 


0.8535343888 


5.0 


1.821763620 


1.28400188 


1.2839939077 


1.2840020996 


10. 


3.199301690 


5.12935962 


5.1293596329 


5.1293596444 


15. 


4.511297510 


12.4337232 


12.433723259 


12.4337232589 


20. 


5.805158110 


23.1467952 


23.146795143 


23.1467951431 


30. 


8.359177000 


54.1918175 


54.191817437 


54.1918174372 


40. 


10.8899708000 


97.83692290 


97.836923003 


97.8369230031 



As suggested by Table[8]the series behaves well for large p, it is found that A(R,p) yields a stable series 



in powers of 1/p. To within 0(l/p 7 ), the expansion for (BBT i is: 

A(R, P ) = ( P+ 1)^2B-M±i> + < 2R+1 > (16^ + 40^ + 23) 



+ 



+ 



+ 



4p 4p 2 64p 3 

(32 i? 2 + 68 R + 41) (64 R 3 + 576 # 2 + 1108 R + 681) 

64^ 512p 5 
(256 i? 3 + 1432 R 2 + 2566 i? + 1593) 
512p 6 

(1280 R A + 28160 R 3 + 123680 R 2 + 210448 R + 131707) 

16384 p 7 

(8192 R A + 95040 R 3 + 358368 R 2 + 587512 R + 371061) 

16384 p 8 

(7168 R 5 + 313600 R A + 2607232 i? 3 + 8854496 R 2 + 14149364 + 9039151) 

131072 p 9 

(65536 R 5 + 1366016 R A + 9200576 R 3 + 29011472 i? 2 45621790 R + 29559559) 



(38) 



131072p 



10 



The coefficients up to 0(p 6 ) have been previously published [40] but our computer algebra programs 
allow us to go much further. 



3.6 Other Bases - Algebraic Vindication 

Although our previous results are apparently independent of the size of the chosen basis, we must con- 
sider other bases. For the ?] coordinate, we consider the Baber-Hasse and the Wilson bases [14] which are 
described as follows. 



Baber-Hasse: 

l=m 

The recurrence relation is given by: 

% + 3)^ [M£ + + Ql]ae+1 + ai(k)ae + W^T) (Ql ~ 2p£)ae ~ 1 = (40) 

where for m = 0: 

ai (k) = A - p 2 + £(£+ 1) 
a_! = . 

Wilson : 

M(r?,0) = e ^ e w(l - ?? )™/ 2 ^(-l) fc c k (1 - V ) k (41) 

k=0 

The recurrence relation is: 

2(k + 1) (fe + m + 1) c fc+ i + <Ti(fc) c fc + 2[Qi + p (& + m) c k _ x = (42) 

where 

ai(k) = A - p 2 + (m + 1) (m + 2p) + A; (fc + 2 m + 4p + 1) 
c_i = . 



Both of these bases have been implemented into the Maple system. If we consider £ = 0, the coefficient 
apf of Baber-Hasse basis is of order 0(l/p ) and the coefficient cn of Wilson basis is of order O(p ). 
However, if we inject our series solution for A(p 2 ) into the series coefficients of both bases, we find that 
both otv and cn are formally zero to within order 0{p N ). This can be seen through a number of computer 
algebra demonstrations. Thus our series solution for A{p 2 ) also formally satisfies the recurrence relations 
of these other bases, order by order in p. 

For the £ coordinate, apart from the used Hylleraas basis, there is also the Jaffe basis. 



Jaffe : 

ho = (e - i) m/2 (e + i)-m-i+Q,/2v e -pt D k (i=±) (43) 

The recurrence relation is: 

(k + 1) (k + m + 1) D k+1 + 7i(fc) D k + {k-^){k + m- ^)D k _ x (44) 

2p 2p 



where 

71 (fc) = A - p 2 + Q 2 - (m + 1) (2p + l 



2p 



2k (k + m + 2p + l- —) 
2p 



= . 

Similarly, it can be algebraically demonstrated that e.g. for j = n — t — 1 = 0, the 1/p series expansion 
of A(R, p) formally satisfies the coefficients of the Jaffe basis for negative powers of p just as they satisfy 
the Hylleraas basis. This demonstration allows us to consider another basis of importance for the r\ coor- 
dinate, namely the Power basis: 



Power : 

M(r], 4>) = e im ^e~ q ^ 1+ri \l — r\) m l 2 d k M(-(k + 8k), m + 1, 2p(l + rf)) (45) 

fc=0 

M is the confluent hyergeometric function. The recurrence relation is: 

(k + 5k + l)(k + 5k + l-% 4+i + Xi (*0 d k (46) 

2p 

Qi 

+ (k + 5k + m)(k + 5k + m ) d k -.\ 

2p 

Xi (fc) = A - p 2 - Q x + (m + 1) (2p- 1 + 

2p 

0i 

- 2(fe + <5A;)(A: + (5fe + m + l-2p-— ) 

2p 

d-i = 0, 

and <5fc is an exponentially vanishing term in R and consequently we do not make the same demonstration 
as for the Wilson and Baber-Hasse bases. However, when we let R — ► oo then 5k — ► and we can make 
a similar demonstration as for the Jaffe basis using the 1/p expansion of A(R,p). 



where 



Granted, we have not proven this for all bases. Nonetheless, we emphasize that e.g. the Wilson basis 
is very different from the Baber-Hasse basis or the Power Basis and the basis of spherical harmonics we 
used as a starting point for this analysis. Moreover, the Hylleraas basis is also very different from the Jaffe 
basis. These demonstrations strongly suggest basis independent results for A{p 2 ) and A(R,p). 

This analysis herein exploits the fundamental theorem of algebra i.e. that if one knows all the N roots 
of a given polynomial say Pjy(x), the latter is completely defined within a scaling factor namely the 
coefficient of its highest power in x. The three-term recurrence relations of eqs. dT9t . d20t and (fT8t have 
a linear dependence on A for the term in dk but no dependence of A for the third term in Thus, 
det(J r ej), det(J r Oj) and det(J^) are i degree polynomials in A regardless of whether or not the third 
term in dk-i is neglected. This allows us to completely account and identify the the eigenparameters of 
the matrices T and y for every discrete state. 

3.7 Mathematical Classification of Solutions 

So far, we have identified the functions implied by the determinants det(J r j) and det(3 ; i)' namely A(p 2 ) 
and A(R, p) respectively for all discrete states where m = for the homonuclear case. In view of previous 
and recent work on the D — > 1 version of and the findings in this work concerning the D = 3 version 
of Hj, we are now equipped with the means to make the following analytical comparison. Here, we can 
put the D — ► 1 and the D = 3 versions of on the same "canonical" footing: 

D — > 1: To reiterate the results of section|2j the energy eigenvalue is governed by an equation of the form: 

d 2 

exp(-2Rd) = P 2 (d) {PN(d)} where E = -— (47) 

When the second order polynomial P2(d) factors into a product of first order polynomials, both 
sides of eq. d47l factors and the solution for d is a (standard) Lambert W function [19,20]. When it 
does not factor, the solution is a generalization of the W function reported in the work of [28]. When 
the right side is a polynomial, the solution is a generalized Lambert W function [30]. The subscript 
Pn (d) reminds us that our generalization for the W function can accommodate a polynomial with 
rational coefficients of arbitrary degree on the right side of eq. d47b . 

The exponential term on the left side is a reflection of the fact that outside the Dirac delta function 
wells, the basis of the particle is a combination of free particle solutions which required matching 
at the Dirac delta function peaks. 

D=3: To summarize the results of the past few sections, the eigenparameter p, which plays an analogous 
role to the parameter d of the D — ► 1 version of H^" is determined from the equation: 

2 

A(R,p) = A( P 2 ) {Pn[p)} where E = -2 (48) 

The subscript Pn{p) reminds us that we have a Taylor series for A(p 2 ) with rational coefficients 
which exactly matches the generalized right-hand side form of eq. (147 1 . However, the left-side 
of d48t looks very different than the left side of (147b : it is the function implied by det(J^) and 
is associated with the separated quantum number j = n — I — 1. Nonetheless, like exp(— 2Rd), 
this function is well-defined asymptotically for large R. The right side of eq. (148 1 is implied from 
det(J r j) for even or odd I which is a united atom quantum number. 

So far, the functions A(R,p) and A{p 2 ) appear in the literature as expansions in terms of p~ k and 
p k respectively, restricted to k = 6 and for specific cases of large and small values of R [38-40]. We 
can obtain series representations of both to a much greater extent in view of our computer algebra 



implementations. We have also seen that A(p 2 ) can be represented as an infinite series in x where 
x = p 2 and is consequently polynomial-like. 

Note that if m ^ 0, the governing equation has the same form as eq. d48l but the left side is more 
complicated and more difficult to get, as the determinant det(J^i) is no longer governed by a simple 
recurrence relation. However, in principle, eq. d48l governs the entire homonuclear case. 

Mathematically, in both cases, the right side of the governing equation is expressed in terms of only one of 
the eigenparameters whereas the left side requires the parameter and the value of the internuclear distance 
R which is determined on input. 

We therefore come to the conclusion that the eigenparameter p, like its D — ► 1 counterpart d, is also deter- 
mined by a special function which is an implicit function, an even greater generalization of the Lambert W 
function. So far, the functions A(R,p) and A{p 2 ) do not appear in the literature. However, we can obtain 
series representations of both to the extent of getting reliable numbers, as demonstrated by our tables of 
values. 

On the subject of implicit functions or implicit equations, these are seen in a number of specific contexts: 

Retardation Effects: Equations of form e~ c x = ^(A) and more generally e~ cX = P2{\)/Q\{X) ex- 
press the solutions of a huge class of delayed differential equations [46, eq.(3)]. 

Bondi's K-calculus: It is well known in the area of special relativity that the Lorentz transformation can 
be derived from the theory of implicit functions with minimal assumptions of continuity [47]. Here 
one seeks the function fit) satisfying f(f(t)) = k 2 t and the requirements that it be monotone 
increasing and continuous. The unique solution is: 

f(f(t))=k 2 t where v/c = k 2 - 1/k 2 -> f(t) = kt 

GRT/QFT: As we mentioned before, the Lambert W function and its generalization appear in General 
Relativity as solutions to respectively the two-body and three-body linear (1 + 1) gravity problems 
via dilaton theory [28, 29]. 

Implicit functions often appear in problems with retardation effects, relativistic or otherwise. Thus, with 
some reservations, we associate with this mathematical category a tentative "physical picture": 

Although the hydrogen molecular ion H^ in the context of the Schrodinger wave equation is not a rela- 
tivistic formulation, the eigensolutions we obtain nonetheless suggests something akin to a retardation or 
delay effect. This is not the case for a one center problem like the hydrogen atom but this characteristic 
appears for a two-center problem. However, this statement must be tempered with the fact that e.g. the 
Lambert W function also appears in many other types of problems with no relationship to retardation 
effects. 



4 Summary/Conclusions 

Through experimental Mathematics using computer algebra as a tool, we have identified the mathemati- 
cal structures governing the energy spectrum of the hydrogen molecular ion H^~ for the two-center one- 
electron problem. 

In the present work, we started with a particular choice of basis and expressed the determinantal conditions 
by which the eigenparameters p and A are obtained. From one of the two determinants, we inverted the 



problem to obtain a series representation of the separation constant A(p 2 ) associated with the united atom 
quantum number I. We applied a similar approach to obtain A(R,p) from the remaining determinant and 
associated with the separated atom quantum number j = n — I — 1. where n is a united quantum number. 
We then demonstrated that the results were independent of the size and even the choice of basis. 

The eigenparameter p for which E = —2p 2 /R 2 is obtained by matching A{p 2 ) = A(R,p) and found 
to be the solution of an implicit function, with features similar to that of the Lambert W function and its 
recent generalizations [30]. This allowed us to mathematically categorize the eigenvalues (or rather make 
us realize what they are not) and even to associate a tentative "physical picture" to the solutions. While 
we made no pretense at rigor, the solutions were nonetheless vindicated numerically and by algebraic 
demonstrations with computer algebra. 

The results express analytical solutions for the ground state and the countable infinity of discrete states of 
for the homonuclear case when the magnetic quantum number m = 0. From the discussion below 
eq. d48b . we anticipate that the eigensolutions for m / for the homonuclear case to be qualitatively 
similar though admittedly this remains to be proven. We emphasize that although the basis and approach 
used here were ideal for m = and the homonuclear case, the computer algebra methods shown are 
directly applicated to the heteronuclear case with m / 0. For m ^ 0, one should work directly from the 
recurrence relations of the chosen basis now that we understand how these basis coefficients behave with 
better and better accuracy for the series expansions of A(p 2 ) in p and the asymptotic series expansions for 
A(R,p) in 1/p. 

However, we make no pronouncements concerning the nature or mathematical category of the solutions 
for the heteronuclear case or when the nuclei are allowed to move. We note that for m = 0, the matrix 
for det(3^) remains tridiagonal while the band matrix for det(^ r j) is pentadiagonal and consequently 
governed by nested recurrence relations [15] suggesting that the analysis shown herein is possible. 

A number of issues arise from this result. In a sense, the result is both overdue and premature. It is 
overdue because of our present capacity to find solutions to fair-sized molecules using computational 
chemistry. On the other hand, it is premature. The functions we found A(p 2 ) and A(R,p) do not seem 
to resemble anything we have seen in the literature. The apparent singularities or "resonances" at 2p 2 = 
3 (4 j + 3) (4 j + 7) for A(p 2 ) and R = (j + 2) p + 2 p 2 for A(R, p) for j = 0, 1, 2, . . . do not constitute 
a problem since the eigenparameters p and A for a given R are never found on these resonances. Once 
a value of R is injected into A(p 2 ) = A(R,p), solving for p numerically did not create any problems 
in the test cases examined so far. At any rate, the tables shown herein merely illustrate the convergence 
properties of the functions A{p 2 ) and A(R,p) we have identified: solving the coupled set of polynomials 
det( Ti) and det(^j) for p and A at a given distance R involves no resonances and is still the most useful 
method from a computational point of view. In principle, the latter can go further than any FORTRAN 
program. 

We have ordered series representations to relatively high order of both of these functions A{p 2 ) and 
A(R,p) and we can generate reliable numbers for a number of discrete quantum states. We have also 
demonstrated that we could use these series beyond their radius of convergence using techniques for han- 
dling divergent series. 

From here, one could explore and seek alternate representations of these functions with better convergence 
properties especially at low R for A(R,p) and large R for A(p 2 ) but the results from the Sidi transfor- 
mations are already very promising. At any rate, the hydrogen molecular ion for clamped nuclei can be 
entirely contained within simple computer algebra sessions, not much more complicated than those of the 
hydrogen atom 1 . 

The exploratory and roundabout way by which we found our solutions, suggests there is something miss- 
1 Maple CAS programs are available upon request. 



ing in the mathematical physics or the methods for obtaining the eigenvalues of the Schrodinger wave 
equation. There is hardly any existing "technology" for solving quantum chemistry problems involving 
implicit functions. Our use of a basis is certainly valid to demonstrate or prove a result. Furthermore, 
the convergence of the bases used here has been confirmed by determining the asymptotic behavior of 
the expansion coefficients of the wavef unctions for the various basis sets considered [40]. Nonetheless, a 
more direct way of generating the functions of A(p 2 ) and A(R,p) would be instructive. 
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Energy vs. R for (a.u.) 

D -> oo D^l D = 3 
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E 3 (R) « ^Ei (f ) + (^) 
Reference: Lopez-Cabrera, Tan and Loeser, 7. Pfryj. C/iem. 97, 2467-2478 (1993). 



